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Abstract. Using an adiabatic approximation we derive an effective interaction 
potential for spatially indirect excitons in quantum well structures. Using 
this potential and path integral Monte Carlo simulations we study exciton 
crystallization and the quantum melting phase transition in a macroscopic system 
of 2D excitons. Furthermore, the superfluid fraction is calculated as a function 
of density and shown to vanish upon crystallization. We show that our results 
complement the previous studies of quantum dipole systems: We predict a 
second quantum phase transition - melting of the crystal at high densities which 
originates from softening of the short-range part of the inter-exciton interaction. 
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1. Introduction 

In recent decades systems of indirect excitons have been extensively studied 
experimentally [1] due to the prospects to achieve superfluidity and Bose condensation. 
Promising setups which have been successful in controlling the many-exciton state are 
electron-hole bilayers, e.g. [2 [3] or single quantum wells using the quantum Stark 
confinement (QSC). The latter will be studied in the present work. As was shown in 
Ref. [3] trough the QSC one can simultaneously produce spatially indirect excitons and 
achieve their spatial localization. By varying the applied electric field one can control 
the exciton coupling parameter - the ratio of the mutual exciton-exciton interaction 
to the confinement energy, as shown in our previous simulations for GaAs and ZnSe 
structures [4] . With the increasing of inter-exciton coupling a striking phenomenon is 
expected - spatial ordering of excitons into a crystalline lattice. This quantum phase 
transition has been recently studied by a path integral Monte Carlo (PIMC) method 
for trapped finite Coulomb systems [5] , systems of particles with dipole interaction [6] 
and for the symmetric electron-hole bilayers [71 [21 [51 [3] . 

In the case of interacting Bose particles, in addition, (partial) superfluidity is 
expected. Then during the crystallization transition one expects a gradual decrease of 
the superfluid fraction. This behaviour has recently been studied in detail by PIMC 
simulations for the Bose particles with the Coulomb interaction [9]. It was shown 
that the superfluid density can be concentrated either in the core or at the cluster 
boundary depending on the hexagonal order in the clusters which sensitively depends 
on the particle number. Furthermore, a two-dimensional homogeneous dipole system 
of up to 400 particles has recently been studied [TOl [11] . The authors found a melting 
point at a density of nr^ — 290 ± 30 or D = \/n'^o = 17 ± 1, for the definitions see 
Eq. ([T4|) . and a step- like vanishing of the superfluid fraction at this point. In the 
present paper we test the validity of the dipole model for indirect excitons. In fact the 
effective exciton-exciton interaction Vxx is a quantum-mechanical four-body problem 
which has been studied e.g. by Zimmermann [T7 and may significantly deviate from 
a simple dipole form at distances comparable with the exciton dipole length. We 
therefore, devote special care to compute Vxx from first principle PIMC simulations. 
We find that this interaction allows for a second phase transition - quantum melting of 
the exciton crystal by compression - missing in dipole systems. Further, we found that 
for an effective exciton dipole length below d < 6a^ the crystalline phase vanishes. 

Below we present results for a model similar to the experimental setup of the 
Timofeev group [13]: a single ZnSe-based QW with indirect excitons produced by 
an electric field [1] applied normal to a QW plane. In order to be able to perform 
simulations for a macroscopic ensemble of indirect excitons we apply a bosonic model. 
As shown by various groups [14[ 115] . in the moderate density regime the excitons 
are adequately treated as a composite particles obeying Bose statistics due to strong 
attractive interaction between electrons and holes. In the simulations of Ref. [15] we 
were able to test this approximation against the exact two-component fermion system 
in a broad range of densities. In particular, we found that the bosonic model gives 
accurate predictions for the superfluid fraction once the excitons are in a strongly 
coupled - low/moderate density regime. In contrast, in the weakly coupled - high 
density regime close to or beyond the Mott density, the results are strongly affected 
by the Fermi statistics, and the bosonic approximation is no more valid. Thus, in the 
present analysis of moderate-density systems the bosonic model works well allowing 
us to study relatively large exciton ensembles without being hampered by the fermion 
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sign problem. 



2. Model 



The general Hamiltonian for the system of Ne electrons and Nh holes [N — 
in the quantum well confinement and E-field can be written as 

with the single particle and interaction parts defined as 
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We consider a homogeneous electric field Ez{ri,Zi) = Ez{zi), vectors r denote 2D 
vectors in the QW plane, V*^^ is the QW confinement, is the electrostatic 
potential energy due to the electric field and e is the background dielectric constant. 

To apply the approximation of bosonic excitons valid for low to moderate 
densities ^5j, we want to reduce the 3D Hamiltonian ([I]) to a 2D one, where all 
eflFects related to a particular width of the QW and the electric field strength will 
be combined in an effective inter-exciton interaction Vxx{R)- This becomes possible 
by using the adiabatic approximation [16|, I17j . This approach is justified for high 
values of the electric field Ez which leads to a strong localization of electrons and 
holes at opposite edges of the QW. To be specific, the calculations below correspond 
to i — 30 . . . 120nm wide ZnSe QWs and E^ — 20kV. Assuming the relation of 
energy scales, AeJ™^'"^ ^Ihlee hh)^ where Ae^*"^'*^ is the characteristic spacing of the 
quantized one-particle energy levels in z-direction und the interaction energy, we 
separate the out-of-plane motion and solve the 3D Bloch equation for the A^-particle 
density matrix p^o 
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in the adiabatic approximation, i.e. 



P3d{P) = P2Diri,---,rN,P) Y[peizt,/3) Y[ph{zj,(3), 



(3) 



(4) 



where /? = l/fc^T (below we drop the argument /?). Now integrating out in Eq. ([3]) all 

z-dependencies, i.e. applying J Y[ dzi Y[ dzj, we obtain a reduced 2D Bloch equation 
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where we introduced a smoothened Coulomb potential 

(r-y ) = j Pa. {zi) {zj) dz, dzj, a,P = e,h. (6) 

The densities Pe(-Ze) and ph{zh) are found by solving a single-exciton problem in an 
electric field |T6]. Also the exciton dipole moment, ^ = e ■ d, follows directly from the 
electron and hole densities 

d^{z.)-iz,)=Jz.p.iz.)dz.-Jz,p,iz,)dz,. (7) 

The considered here low to moderate density regime leads also to another relation 
of energy scales, i.e. Eb{X) ^ YxxiksT, where Eb is the exciton binding energy. 
Under these conditions the excitons remain in their internal quantum states described 
by a two-body density matrix p'^^{Yeh) throughout their interaction. This pair density 
matrix depends on the electron-hole separation r^h = J^e — 'Th and can be obtained 
numerically, e.g. with the matrix-squaring technique |18| applied to the interaction 
potential Using again the adiabatic approximation (now in the 2D plane) we write 
p2D as a product of a density matrix of an particle complex and relative density 
matrices p^^ of excitons 

P2D = P2d(R\...,R^^) Up'^t:,). (8) 

a=l 

Here we have assumed electrical neutrality, Ne — Nh — Nx, and introduced the 
electron- hole pair coordinates related to the same exciton, (ri, . . . ,rjv) = {{r'^,rf^) = 
(R°, r°^)}|o=i I with the center of mass (c.o.m.) coordinates R" = (mUrg -|- 
m|[r^)/Ma;, Mx = ml + mj[. Certainly, the ansatz ([8]) implies that the excitons are 
stable against external perturbations and we are below the Mott density. 

Now averaging Eq. ([5]) over the relative degrees of freedom of excitons, i.e. 

integrating over J Y[ dv^^, we obtain the 7Vj;-exciton Bloch equation depending on 

a=l 

the c.o.m. coordinates 

9p2d(R\...,R^^) 



a/3 
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The interaction term in the effective exciton Hamiltonian H'^^^ is defined as the sum 
of the effective (adiabatically averaged) interactions of two electrons and two holes in 
excitons o and b {a ^ b) 

Vxx{R''')^j Vo.p{K-A\)P"''i<h)p'''i<h)d<hdvlH- (12) 

Q,/3— e,/i 

The distances of two particles from different excitons can be expressed as 
r^, - r^, R'^ - R^ + m^iv^,, - r^J/M,, 

- = R'^ - R" - mu{vl^ - r^J/M„ 

- = - R" + (m.r^;, -f m;.r^J/M,. (13) 
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Table 1. Semiconductor QW parameters. Masses are in units of the free electron 
mass mo. 





GaAs/AlGaAs 


ZnSo/ZnSSo 


e 


12.58 


8.7 


m| 


0.0667 


0.15 




0.112 


0.37 




0.377 


0.86 




2.68 


3.46 




9.98 


3.07 


Ha* [mcV] 


11.47 


53.93 



After integration, in Eq. (fT^ remains only the c.o.m. dependence on = jR" — R''|. 

Thus, we have derived the effective Hamiftonian pil)l of composite particles. The 
interaction potential (fT2|l generalizes the dipole potential used in the previous analysis 
of spatially indirect excitons [6l [101 E] ■ The comparison of both is discussed below. 
The corresponding iV-body problem ([T0|) can be solved with the path integral Monte 
Carlo technique which allows for a direct treatment of many-body correlation and 
bosonic exchange effects, for details see [111 HOI HI] ■ 

3. Results 

We have performed PIMC simulations for a 2D homogeneous system with iV = 60 
and 90 bosonic excitons in a simulation box with periodic boundary conditions. The 
potential (fT^ has been divided into a short and long-range part, V^x = ^xx—yD)+yD , 
with the dipole interaction, Vd = {ed)'^/er^, treated by the usual Ewald summation 
technique. We used the following system of units: r ^ r/a*g, E ^ E/Ha*, with the 

electron Bohr radius, Og = ft^e/meC^, and the electron Hartree, Ha* — e'^/ea%. Here 
Mx = ml + mj[ is the exciton mass. Parameters for typical semiconductor structures 
are listed in table [H As it follows from the derivation in adiabatic approximation the 
Hamiltonian (|10p contains only the in-plane particles masses, Tijjj^j. The anisotropy of 
the parabolic bands both for electrons and holes, i.e. the out-of-plane effective masses, 
m^,, N, are involved in the exciton solution in the z direction. These masses determine 
the shape of the density matrices pe{ze), Phizh) and hence indirectly influence the 
effective inter-exciton interaction via Eq. ([6]). In particular, for a 20kV/cm electric 
field applied to a 30nm wide QW the calculations of Refs. [H [TB] predict that the e-h 
separation ([7]) equals d = 15.78nm for GaAs and d = 20.41nm for ZnSe structures. 
Now comparing the ratio d = d/a*g for both structures, we find d = 1.58 and d = 6.65, 
respectively. This shows that in ZnSe QW the excitons are more strongly coupled and 
it is easy to reach a crystalline regime as discussed below. Other advantages of using 
materials with larger effective masses are: a) increased stability of excitons due to 
higher binding energies and b) increased exciton life-time due to a better separation 
of carriers in the z-direction [the radiative life-time depends on the overlap of the 
density matrices Pe{h){ze(h))]- 

In the following we, therefore, concentrate on the ZnSe structure. However, the 
results presented in the dimensionless units, r/a*g and E/Ha*, using table [1] can be 
applied to other materials as well. In particular, the effective potential Vxx (fT2|) already 
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Figure 1. (Color online) Left: exciton interaction potential Vxx{'i')[Ha*\, 
Eq. 1 II2I I. for several dipole moments d. Right: Vxx/d? compared with the 
dipole potential, VuiT) = l/f^, and classical exciton potential, Vea; = (2/?* — 
2/\/r'^ + d?)/d^ (shown for d = do)- 




reduces to a dipole interaction Vd = (ed)'^ /er^ at distances of about several exciton 
dipole moments d. In this case the Hamiltonian (|10[) can be brought to a universal 
dimensionless form using the scale [11]: oq = Xj^/n. Eq ~ fi^ /MxO^. Both systems of 
units are connected via relations 

d — d/a*g^ 

(14) 

where n — N/{LxLy) is the number density. 

In figure [1] we show Vxx{''')[Ha*] for several e-h separations, d — do,2dQ,3dQ,AdQ 
with do = 6.64848 [d — do corresponds to a 30nm ZnSe QW and field strength 
20kV/cm]. In the right part we show, in addition, the dipole potential Vd and the 

classical exciton pair potential (limit of d ^ r). Vex = 2/r — 2/ + cP. While at 
large distances Vxx agrees with Vd, for r < 4c?, Vxx is substantially weaker. Further, 
as one might expect from the dipole model, overall the interaction is stronger with 
increasing d. However, at small distances, r < 6a|j, Vxx shows the opposite behavior 
which originates from the smoothening procedure (fT2|) over the exciton relative density 
matrices: with increase of d e-h pairs become more weakly bound, and the exciton 
in-plane size increases. This delocalization reduces the strength of the Coulomb 
interaction between two electrons and two holes, i.e. between two excitons. Vice 
versa, the stronger the binding of an e-h pair and its spatial localization the faster Vxx 
approaches Vex- 

Let us now analyze the melting behavior of the exciton ensemble. In Refs. pUilllj . 
crystallization of dipoles was observed at D = 17 ± 1 which was identified by a 
break of the translational symmetry of the pair distribution functions (PDF) and the 
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static structure factor. We now perform a similar analysis with the improved model 
interaction Vxxii")! for T = l/3000Ha, and consider a density range of pa*^ — l/nr'^, 
rs = 5 . . . 12. Simulations were performed for d[do] = 1,2,3,4, and the results 
are shown in figures [5] and [31 At low densities where our potential is close to a 
dipole potential we observe similar results as in [101 111] , i.e. crystallization of dipoles 
upon compression (not shown). But most importantly, at high density, we observe 
completely different behavior which is due to the weak potential at small r: the 
exciton crystal melts upon compression. This transition is clearly seen from the 2D 




Figure 2. (Color) Left: 2D pair distribution function at rs = 10.0 (top) in the 
solid and rs = 9.5 (bottom) in the supcrfluid gas phase. Right: snapshots from 
PIMC simulations. Trajectories of particles involved in permutations of different 
length are denoted with different colors, d = 3do 

PDF in figure [2l^left) where an abrupt loss of a (quasi) long-range crystalline order is 
observed by a slight change in r^ from 10 (top) to 9.5 (bottom). Simultaneously, in 
figure[Hright) we observe a topological change in the picture of the particle trajectories 
in the path integral representation [19]. While in the solid phase we observe only 
local exchanges of few particles, just after the melting transition the trajectories form 
macroscopically large permutation cycles crossing the edges of our periodic simulation 
cell. From the statistics of the flux of paths winding around the periodic cell one can 
estimate the fraction of the superfluid density [THl [H] : 

7, = p,/p = Mx{W^)lfepNx, W = Y, I dt [dv,{t)/dt\ . (15) 

Figure [3] indicates a step-like increase of the superfluid density from zero up to about 
35% in the gas phase. The vertical dotted line shows the Mott density where the 
excitons pressure ionize and the bosonic model fails. The critical values for the 
exciton quantum melting transition at different d are collected in table [2] together 
with the critical data r\ for the dipole crystallization. 

In conclusion, the derived exciton-exciton potential leads to completely different 
predictions for the phase diagram of bosonic excitons compared to the dipole model. 
Due to the much softer Coulomb-like interaction at small distances, the exciton solid 
melts by compression, similar to a Wigner crystal of electrons ^ . Due to this fact it 
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Table 2. Interparticle distances at first and second phase transition; 
superfluid gas — exciton solid. r\ are estimated from II14I I and D > 17 ± 1; are 
the PIMC results using Vxx{t) in figure [T] The exciton solid exist for densities 
na*^ = l/vrr^, with < < '"i- Second colum, L, is the required ZnSe QW 
width at the field strength = 20kV/cm . 



d[do\ 


L [nm] 


' S 


' s 


1 


30 


5.1(0.3) 




2 


~ 50 


20.4 (1.2) 


10.0 (0.5) 


3 


~ 70 


45.90 (2.7) 


10.0 (0.5) 


4 


- 90 


81.6 (4.8) 


11.0 (0.5) 



becomes possible to stabilize the exciton lattice only in a finite density interval [see 
table [2]. Outside of this region the excitons exist in a superfluid gas phase. 

Several heterostructures are candidates for the observed effect, but ZnSe is 
favorable due to its relatively high value of the dipole moment. Using parameters 
from table [1] we estimate the exciton solid to exist in ZnSe (taking d — 2(io) in a QW 
with L - 50nm between 0.81 < p[10^°cm-'^] < 3.38 at T <2K and in GaAs in a QW 
with L - 148nm between 0.77 < p[10°cm-'^] < 3.2 and T < 0.4K. While we have 
not considered disorder effects due to the imperfections of the QW planes they can 
be important. In our case of the electric field- induced indirect excitons electrons and 
holes are pushed to the QW edges and hence experience the influence of QW width 
fluctuations and impurities [17] . As in the case of the electron Wigner crystal this can 
additionally stabilize the exciton solid at high densities. 



0,25 




0,000 0,005 0,010 

Density, na^-(n r^) ' 

Figure 3. Supcrfiuid fraction vs. density na*^ , for Af = 60 (circles) and = 90 
(rombs) for d = 3do- Dotted line shows the Lindemann parameter Ur- Vertical 
line indicates the Mott density. The Bose model (superfiuid data) is applicable 
only at lower densities |15| . 
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